class: center, middle, inverse, title-slide # A review of statistical computing with TMB ### Andrea Havron
NOAA Fisheries, OST --- layout: true .footnote[U.S. Department of Commerce | National Oceanic and Atmospheric Administration | National Marine Fisheries Service] <style type="text/css"> code.cpp{ font-size: 14px; } code.r{ font-size: 14px; } </style> --- # Statistical Computing Review <br> .pull-left[ <img src="data:image/png;base64,#static/TMB-trifecta.png" width="60%" style="display: block; margin: auto auto auto 0;" /> ] --- class: middle ### Review: --- # ML Inference #### What is the likelihood of getting 30 heads in 100 coin flips? .pull-left[ 1. **Specify the model** <br><br> `\(y ~ \sim Binomial(n, p)\)` ] --- # ML Inference #### What is the likelihood of getting 30 heads in 100 coin flips? .pull-left[ 1. Specify the model 2. **Calculate the likelihood**<br><br> `\(L(p; n, y) = \frac{n!}{y!(n-y)!}p^y(1-p)^{n-y}\)` ] .pull-right[ `\(y = dbinom(x = 30, n = 100, p)\)`
] --- # ML Inference **What is the likelihood of getting 30 heads in 100 coin flips?** .pull-left[ 1. Specify the model 2. Calculate the likelihood 3. **Calculate the negative log-likelihood**<br><br> `\(-\ell(p; n, y) = -[ln\big(\frac{n!}{y!(n-y)!}\big) + yln(p)\)`<br> `$$+ (n-y)ln(1-p)]$$` ] .pull-right[ `\(nll = -log(dbinom(x = 30, n = 100, p))\)`
] --- # ML Inference What is the likelihood of getting 30 heads in 100 coin flips? .pull-left[ 1. Specify the model 2. Calculate the likelihood 3. Calculate the negative log-likelihood `\(-\ell(p; n, y) = -[ln\big(\frac{n!}{y!(n-y)!}\big) + yln(p)\)`<br> `$$+ (n-y)ln(1-p)]$$` 4. **Calculate the derivative w.r.t. `\(p\)`**<br><br> `\(\frac{d(\ell(p; n, y))}{dp} = \frac{y}{p}- \frac{n-y}{1-p}\)` ] .pull-right[ `\(-ln\big[L(p; n = 100, y = 30)\big]\)` <img src="data:image/png;base64,#00_00_Statistical_Computing_Review_files/figure-html/unnamed-chunk-4-1.png" width="80%" /> ] --- # ML Inference What is the likelihood of getting 30 heads in 100 coin flips? .pull-left[ 1. Specify the model 2. Calculate the likelihood 3. Calculate the negative log-likelihood 4. Calculate the derivate wrt p 5. **Set to 0 and solve for MLE**<br> `\(0 = \frac{y}{p}- \frac{n-y}{1-p}\)` <br> `\(E[p] = \frac{y}{n}\)`<br> `\(E[y] = np\)` ] .pull-right[ `\(-ln\big[L(p; N = 100, y = 30)\big]\)` <img src="data:image/png;base64,#00_00_Statistical_Computing_Review_files/figure-html/unnamed-chunk-5-1.png" width="70%" /> <br> `\(\hat{p} = \frac{30}{100} = 0.3\)` ] --- # ML Inference What is the likelihood of getting 30 heads in 100 coin flips? .pull-left[ 1. Specify the model 2. Calculate the likelihood 3. Calculate the negative log-likelihood 4. Calculate the derivate wrt p 5. Set to 0 and solve for MLE<br> `\(0 = \frac{y}{p}- \frac{n-y}{1-p}\)` <br> `\(E[y] = np\)` 6. **Approximate Var[p] using the second derivative**<br> `\(-\frac{y}{p^2} - \frac{(n-y)}{(1-p)^2}\)`<br> `\(-\frac{np}{p^2} - \frac{(n-np)}{(1-p)^2}\)`<br> `\(-\frac{n}{p} - \frac{n(1-p)}{1-p}\)`<br> `\(l''(p) = -\frac{n(1-p)}{p}\)`<br> `\(Var[p] = \frac{p}{n(1-p)}\)` ] .pull-right[ `\(Var[p] \approx -\frac{1}{l''(p)}\)`<br> `\(SE[p] \approx \sqrt{ \frac{.3}{100(.7)}} \approx 0.065\)`<br> <img src="data:image/png;base64,#00_00_Statistical_Computing_Review_files/figure-html/unnamed-chunk-6-1.png" width="70%" /> ] --- #Multivariate asymptotics * For N-d models, the curvature is represented by a NxN **Hessian** matrix of 2nd partial derivatives * Inverting the negative Hessian gives us a covariance matrix .pull-left[ `\begin{align} (\mathbb{H}_{f})_{i,j} &= \frac{\partial^2f}{\partial \theta_{i}, \partial x\theta_{j}} = \frac{-1}{Var(\Theta)} \end{align}` <br> [**What causes a singular covariance matrix?**](https://andrea-havron.shinyapps.io/mvnplot/) ] .pull-right[ <!-- --> ] --- # What is AD? .three-column-left[ **Automatic Differentiation**<br> Derivatives calculated automatically using the chain rule<br> .p[ - Efficient: forward mode, O(n); reverse mode, O(m) - Accurate - Higher order derivatives: easy ]] .three-column-left[ **Symbolic Differentiation**<br> Computer program converts function into exact derivative function<br> .p[ - Inefficient: O(2n) for n input variables - Exact - Higher order derivatives: difficult due to complexity ]] .three-column-left[ **Numerical Differentiation**<br> Approximation that relies on finite differences .p[ - Efficient: O(n) - Trade-off between truncation error versus round-off error - Higher order derivatives calculation difficult due to error accumulation ]] --- # Computational Graph (Tape) <br> .three-column-left[ ```cpp //Program v1: x = ? v2: y = ? v3: a = x * y v4: b = sin(y) v5: z = a + b ``` ] .three-column[ <img src="data:image/png;base64,#static/comp-graph.png" width="100%" style="display: block; margin: auto auto auto 0;" /> ] .three-column-left[ ```cpp //Reverse Mode dz = ? da = dz db = dz dx = yda dy = xda + cos(y)db ``` ] --- # Reverse Mode .pull-left[ **Static (TMB: CppAD, TMBad)**<br> The graph is constructed once before execution .p[ - Less flexibility with conditional statements that depend on parameters. - Atomic functions can be used when conditional statements depend on parameters - High portability - Graph optimization possible ]] .pull-right[ **Dynamic (Stan: Stan Math Library, ADMB: AUTODIF)**<br> The graph is defined as forward computation is executed at every iteration .p[ - Flexibility with conditional statements - Optimization routine implemented into executable - Less time for graph optimization ] ] --- class: middle # Type Systems in R and C++ --- # Dynamic vs. Static Typing <br> .pull-left[ **R: Dynamic** .p[ - Type checking occurs at run time - The values and types associated with names can change - Change in type tends to be implicit ]] .pull-right[ **C++: Static** .p[ - Type checking occurs at compile time - The values associated with a given name can be limited to just a few types and may be unchangeable - Change in type tends to be explicit ]] --- # What is Templated C++? <br><br> * Generic programming * Allows developer to write functions and classes that are independent of Type * Templates are expanded at compile time .pull-left[ ```cpp template <class T> T add(T x, T y){ return x + y; } int main(){ int a = 1; int b = 2; double c = 1.1; double d = 2.1; int d = add(a,b); double e = add(c,d); } ``` ] .pull-right[ ```cpp int add(int x, int y){ return x + y; } double add(double x, double y){ return x + y; } ``` ] --- # TMB AD Systems <br> .pull-left[ **CppAD** - [CppAD package](https://coin-or.github.io/CppAD/doc/cppad.htm) ] .pull-right[ **TMBad**<br> - TMBad is available with TMB 1.8.0 and higher ]